# =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
# replicationCode.r
# Code to replicate the results reported in:
# "What to do about Atheoretic Lags"
# Cranmer, Skyler, Douglas Rice, and Randolph Siverson
#
# Last Updated 3.5.15
# =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

# front end matters
rm(list=ls())
library(foreign)
library(MCMCpack)

# source atheoretic lag BMA functions; you'll need to ensure the replication files are in your working directory.
#setwd()
source("albma.R")

# =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
# Courts Example
# =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

data <- read.dta("CourtsData.dta", convert.factors=FALSE)

## Results from Table 1: Baseline litigant signal model.
formula<-CASES~ LIPP + LLEGIS + LCASES + LMQV + ABMCH + warren + dumcivrights + dumcriminal + dumunions + dumecon + dumdueproc + dumjudicial + dumfed + dumenvironment + dumtax
rfj.out <- MCMCregress(formula, b0=0, B0=0.1, burnin = 1000, mcmc = 10000, 
   marginal.likelihood= "Chib95", data=data)
# Plots of Posterior Distributions.
par(mfrow=c(1,2))
plot(density(rfj.out[,2]), main="Policy Priority t-1", xlab="Posterior Distribution", 
   lwd=2, ylim=c(0,.4))
abline(v=0,lty=2)
plot(density(rfj.out[,3]), main="Legislative Attention t-1", 
   xlab="Posterior Distribution", lwd=2, ylim=c(0,.4))
abline(v=0,lty=2)

## Results from Table 2: Distributed Lag litigant signal model.
formula<-CASES~ LIPP + L2IPP + L3IPP + L4IPP + L5IPP + L6IPP + 
    LLEGIS + L2LEGIS + L3LEGIS + L4LEGIS + L5LEGIS + L6LEGIS + 
	LCASES + LMQV + ABMCH + warren + dumcivrights + dumcriminal + 
	dumunions + dumecon + dumdueproc + dumjudicial + dumfed + dumenvironment + dumtax
dlj.out <- MCMCregress(formula, b0=0, B0=0.1, burnin = 1000, mcmc = 10000, 
   marginal.likelihood= "Chib95", data=data)
# Plots of Posterior Distributions for Policy Priority
par(mfrow=c(2,3))
plot(density(dlj.out[,2]), main="Policy Priority, t-1", xlim=c(-6,8), lwd=2, xlab="")
par(new=TRUE)
abline(v=0, lty=2)
plot(density(dlj.out[,3]), main="Policy Priority, t-2", xlim=c(-6,8), lwd=2, xlab="")
par(new=TRUE)
abline(v=0, lty=2)
plot(density(dlj.out[,4]), main="Policy Priority, t-3", xlim=c(-6,8), lwd=2, xlab="")
par(new=TRUE)
abline(v=0, lty=2)
plot(density(dlj.out[,5]), main="Policy Priority, t-4", xlim=c(-6,8), lwd=2, 
   xlab="Posterior Distribution")
par(new=TRUE)
abline(v=0, lty=2)
plot(density(dlj.out[,6]), main="Policy Priority, t-5", xlim=c(-6,8), lwd=2,  
   xlab="Posterior Distribution")
par(new=TRUE)
abline(v=0, lty=2)
plot(density(dlj.out[,7]), main="Policy Priority, t-6", xlim=c(-6,8), lwd=2, 
   xlab="Posterior Distribution")
par(new=TRUE)
abline(v=0, lty=2)
# Plots of Posterior Distributions for Legislative Attention.
par(mfrow=c(2,3))
plot(density(dlj.out[,8]), main="Legislative Attention, t-1", xlim=c(-6,8), lwd=2, xlab="")
par(new=TRUE)
abline(v=0, lty=2)
plot(density(dlj.out[,9]), main="Legislative Attention, t-2", xlim=c(-6,8), lwd=2, xlab="")
par(new=TRUE)
abline(v=0, lty=2)
plot(density(dlj.out[,10]), main="Legislative Attention, t-3", xlim=c(-6,8), lwd=2, xlab="")
par(new=TRUE)
abline(v=0, lty=2)
plot(density(dlj.out[,11]), main="Legislative Attention, t-4", xlim=c(-6,8), lwd=2, 
     xlab="Posterior Distribution")
par(new=TRUE)
abline(v=0, lty=2)
plot(density(dlj.out[,12]), main="Legislative Attention, t-5", xlim=c(-6,8), lwd=2, 
     xlab="Posterior Distribution")
par(new=TRUE)
abline(v=0, lty=2)
plot(density(dlj.out[,13]), main="Legislative Attention, t-6", xlim=c(-6,8), lwd=2, 
     xlab="Posterior Distribution")
par(new=TRUE)
abline(v=0, lty=2)


# Results from Table 3: Atheoretic Lags model
formula <- CASES ~ IPP + LEGIS + LCASES + MQV + ABMCH + warren + dumcivrights + dumcriminal + dumunions + dumecon + dumdueproc + dumjudicial + dumfed + dumenvironment + dumtax

judicial.out<-alreg(formula=formula, 
	data=data, model="ls", 
	method="pooled", lags=c(1,2), 
	by="svalue", time.var=data$year, 
	min.lag=1, max.lag=8,
	ml.type = "Chib95",
	b0 = 0, B0 = 0.1,
	mod.prob.method = "BF")
# Gather posterior means across specifications for Policy Priority.
ipp.triples <- matrix(NA,8,3)
ipp.triples[1,] <- c(summary(rfj.out)$quantiles[2,3], summary(rfj.out)$quantiles[2,1],
   summary(rfj.out)$quantiles[2,5])
ipp.triples[2:7,] <- c(summary(dlj.out)$quantiles[2:7,3], 
   summary(dlj.out)$quantiles[2:7,1], summary(dlj.out)$quantiles[2:7,5])
ipp.triples[8,] <- c(judicial.out$bma[2,1], judicial.out$bma[2,3], 
   judicial.out$bma[2,4])
# Plot the posterior means across specifications.
plot(c(1:8), ipp.triples[,1], ylim=c(min(ipp.triples), max(ipp.triples)), 
   pch=19, xaxt="n", ylab="Posterior Mean", xlab="Measure", main="Policy Priority")
axis(1, at=c(1:8), labels=c("Base", "DL t-1", "DL t-2", "DL t-3", "DL t-4", 
   "DL t-5", "DL t-6", "BMA"))
segments(c(1:8), ipp.triples[,2], c(1:8), ipp.triples[,3], lwd=2)
abline(h=0, lty=2)     
abline(v=1.5)     
abline(v=7.5)
# Gather posterior means across specifications for Legislative Attention.
ipp.triples <- matrix(NA,8,3)
ipp.triples[1,] <- c(summary(rfj.out)$quantiles[3,3], summary(rfj.out)$quantiles[3,1], 
     summary(rfj.out)$quantiles[3,5])
ipp.triples[2:7,] <- c(summary(dlj.out)$quantiles[8:13,3], 
   summary(dlj.out)$quantiles[8:13,1], summary(dlj.out)$quantiles[8:13,5])
ipp.triples[8,] <- c(judicial.out$bma[3,1], judicial.out$bma[3,3], 
   judicial.out$bma[3,4])
# Plot the posterior means across specifications.
plot(c(1:8), ipp.triples[,1], ylim=c(min(ipp.triples), max(ipp.triples)), pch=19, 
   xaxt="n", ylab="Posterior Mean", xlab="Measure", main="Legislative Attention")
axis(1, at=c(1:8), labels=c("Base", "DL t-1", "DL t-2", "DL t-3", "DL t-4", 
   "DL t-5", "DL t-6", "BMA"))
segments(c(1:8), ipp.triples[,2], c(1:8), ipp.triples[,3], lwd=2)
abline(h=0, lty=2)     
abline(v=1.5)     
abline(v=7.5)


## Figure 1: Posterior Model Probabilities for Atheoretic Lags Model and Across Different Max Lag Choices

# Left Panel: Plot posterior model probabilities.
plot(judicial.out$postprob, type="b", pch=4, lwd=2, xlab = "Lag", 
     ylab = "Posterior Model Probability")
	 
# Right Panel: Choose a group of maximum lags.
mxlg <- c(4,6,8,10)
# Set up matrix for output.
plot.output <- matrix(NA, 4, 10)
# Loop through models with different feasible sets.
for (i in 1:length(mxlg)){
  judicial.out <- alreg(formula=formula,
                        data=data, model="ls", method="pooled", lags=c(1,2), 
                        time.var = data$year, by="svalue", min.lag = 1, 
                        max.lag = mxlg[i], ml.type ="Chib95", b0 = 0, B0 = .1,
                        mod.prob.method = "BF")
  plot.output[i,1:mxlg[i]] <- judicial.out$postprob
}
# Plot posterior model probabilities.
plot(1:10, plot.output[1,], type="b", pch=19, col=rgb(0,0,0,alpha=.3), lwd=2, 
   xlab = "", ylab = "", xlim=c(1,10), ylim= c(0,1))
par(new=TRUE)
plot(1:10, plot.output[2,], type="b", pch=19, col=rgb(0,0,0,alpha=.3), lwd=2, 
   xlab = "", ylab = "", xlim=c(1,10), ylim= c(0,1))
par(new=TRUE)
plot(1:10, plot.output[3,], type="b", pch=19, col=rgb(0,0,0,alpha=.3), lwd=2, 
   xlab = "", ylab = "", xlim=c(1,10), ylim= c(0,1))
par(new=TRUE)
plot(1:10, plot.output[4,], type="b", pch=19, col=rgb(0,0,0,alpha=.3), lwd=2, 
   xlab = "Lag", ylab = "Posterior Model Probability", xlim=c(1,10), ylim= c(0,1))

## Figure 2: BMA Posterior Means Across Max Lag Choice
# Choose a group of maximum lags.
mxlg <- c(3,4,5,6,7,8,9,10)
# Set up matrix for output.
plot.output <- matrix(NA, 8, 3)
# Loop through models with different feasible sets.
for (i in 1:length(mxlg)){
  judicial.out <- alreg(formula=formula,
                        data=data, model="ls", method="pooled", lags=c(1,2), 
                        time.var = data$year, by="svalue", min.lag = 1, 
                        max.lag = mxlg[i], ml.type ="Chib95", b0 = 0, B0 = .1,
                        mod.prob.method = "BF")
# Pull BMA estimate on Policy Priority.
  plot.output[i,1] <- judicial.out$bma[2,1]
# Pull lower bound on BMA estimate on Policy Priority.
  plot.output[i,2] <- judicial.out$bma[2,3]
# Pull upper bound on BMA estimate on Policy Priority.
  plot.output[i,3] <- judicial.out$bma[2,4]

}
# Plot average posterior means.
plot(mxlg, plot.output[,1], pch=19, cex=1.5, xlab = "Feasible Set", 
   ylab = "Posterior Mean", xlim=c(min(mxlg),max(mxlg)), ylim= c(0,11))
# Add 95% credible intervals.
segments(mxlg, plot.output[,2], mxlg, plot.output[,3], col="black", lwd=2)
abline(h=0, lty=2, lwd=2, col="gray24")
par(new=TRUE)

## Figure 3: BMA Posterior Model Probabilities across priors
# Identify the set of prior means.
b0 <- c(-5,-2,-1,-.5,-.1, -.01,0,.01,.1,.5, 1, 2, 5)
# Set up a matrix to store output.
output <- matrix(NA,length(b0),6)
# Loop through the models.
for (i in 1:length(b0)){
  judicial.out <- alreg(formula=formula,
                      data=data, model="ls", method="pooled", lags=c(1,2), 
                      time.var = data$year, by="svalue", min.lag = 1, max.lag = 6,
                      ml.type ="Chib95", b0 = b0[i], B0 = .1,
                      mod.prob.method = "BF")
# Extract posterior model probabilities.
  output[i,] <- judicial.out$postprob
}
# Plot posterior model probabilities. Each line represents a particular lag, 
  as indicated by pch.
plot(b0, output[,1], xlim=c(-5,5), ylim=c(0,1), type = "b", pch="1", 
   col=rgb(0,0,0,alpha=.3), lwd=2, ylab="", xlab="")
par(new=TRUE)
plot(b0, output[,2], xlim=c(-5,5), ylim=c(0,1), type = "b", pch="2", 
   col=rgb(0,0,0,alpha=.3), lwd=2, ylab="", xlab="")
par(new=TRUE)
plot(b0, output[,3], xlim=c(-5,5), ylim=c(0,1), type = "b", pch="3", 
   col=rgb(0,0,0,alpha=.3), lwd=2, ylab="", xlab="")
par(new=TRUE)
plot(b0, output[,4], xlim=c(-5,5), ylim=c(0,1), type = "b", pch="4", 
   col=rgb(0,0,0,alpha=.3), lwd=2, ylab="", xlab="")
par(new=TRUE)
plot(b0, output[,5], xlim=c(-5,5), ylim=c(0,1), type = "b", pch="5", 
   col=rgb(0,0,0,alpha=.3), lwd=2, ylab="", xlab="")
par(new=TRUE)
plot(b0, output[,6], xlim=c(-5,5), ylim=c(0,1), type = "b", pch="6", 
   col=rgb(0,0,0,alpha=.3), lwd=2, ylab="Posterior Probability", xlab="Values of b0")

# Identify the set of prior precisions
B0 <- c(.001,.02,.05,.1,.2,.5,.75, 1, 2, 5)
# Set up a matrix to store output.
output <- matrix(NA,length(B0),6)
# Loop through the models.
for (i in 1:length(B0)){
  judicial.out <- alreg(formula=formula,
                        data=data, model="ls", method="pooled", lags=c(1,2), 
                        time.var = data$year, by="svalue", min.lag = 1, max.lag = 6,
                        ml.type ="Chib95", B0 = B0[i], b0 = .1,
                        mod.prob.method = "BF")
  
  output[i,] <- judicial.out$postprob
}
# Plot posterior model probabilities. Each line represents a particular lag, 
  as indicated by pch.
plot(B0, output[,1], xlim=c(0,5), ylim=c(0,1), type = "b", pch="1",  
   col=rgb(0,0,0,alpha=.3), lwd=2, ylab="", xlab="")
par(new=TRUE)
plot(B0, output[,2], xlim=c(0,5), ylim=c(0,1), type = "b", pch="2", 
   col=rgb(0,0,0,alpha=.3), lwd=2, ylab="", xlab="")
par(new=TRUE)
plot(B0, output[,3], xlim=c(0,5), ylim=c(0,1), type = "b", pch="3", 
   col=rgb(0,0,0,alpha=.3), lwd=2, ylab="", xlab="")
par(new=TRUE)
plot(B0, output[,4], xlim=c(0,5), ylim=c(0,1), type = "b", pch="4", 
   col=rgb(0,0,0,alpha=.3), lwd=2, ylab="", xlab="")
par(new=TRUE)
plot(B0, output[,5], xlim=c(0,5), ylim=c(0,1), type = "b", pch="5", 
   col=rgb(0,0,0,alpha=.3), lwd=2, ylab="", xlab="")
par(new=TRUE)
plot(B0, output[,6], xlim=c(0,5), ylim=c(0,1), type = "b", pch="6", 
   col=rgb(0,0,0,alpha=.3), lwd=2, ylab="Posterior Probability", xlab="Values of B0")

# =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
# Modernization Example
# =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
   
rm(list=ls())
data <- read.dta("ModernizationData.dta")
source("albma.R")

## Table 4: Baseline modernization model 
formula <- demd ~ demd_L1 + income_L1 + lpop_L1 + ethfrac_L1 + relfrac_L1 + instab_L1 + oil + ncontig   
# Estimate baseline model.
rfp.out <- MCMCprobit(formula, b0=0, B0=.01, burnin = 1000, mcmc = 10000, 
   marginal.likelihood= "Chib95", data=data)

## Figure 4: ACF Plot
par(mfrow=c(2,2))
plot(acf(data$income, plot=F), cex=1.5, main="Income")
plot(acf(data$lpop, plot=F), cex=1.5, main="Log(Population)")
plot(acf(data$ethfrac, plot=F), cex=1.5, main="Ethnic Fract.")
plot(acf(data$relfrac, plot=F), cex=1.5, main="Religious Fract.")
   
## Table 5: Distributed lag modernization model
formula <- demd ~ income_L1 + income_L2 + income_L3 + income_L4 + income_L5 + income_L6 + income_L7 + income_L8 + income_L9 + income_L10 + 
   lpop_L1 + lpop_L2 + lpop_L3 + lpop_L4 + lpop_L5 + lpop_L6 + lpop_L7 + lpop_L8 + lpop_L9 + lpop_L10 + 
   ethfrac_L1 + ethfrac_L2 + ethfrac_L3 + ethfrac_L4 + ethfrac_L5 + ethfrac_L6 + ethfrac_L7 + ethfrac_L8 + ethfrac_L9 +
   instab_L1 + instab_L2 + instab_L3 + instab_L4 + instab_L5 + instab_L6 + instab_L7 + instab_L8 + instab_L9 + instab_L10 + demd_L1 + oil + ncontig
# Estimate distributed lag model.
dlp.out <- MCMCprobit(formula, b0=0, B0=.01, burnin = 1000, mcmc = 10000, 
   marginal.likelihood= "Chib95", data=data)

# Compose results table.
su.out <- matrix(NA,nrow(summary(dlp.out)$statistics),4)
su.out[,1] <- summary(dlp.out)$statistics[,1]
su.out[,2] <- summary(dlp.out)$statistics[,2]
su.out[,3] <- summary(dlp.out)$quantiles[,1]
su.out[,4] <- summary(dlp.out)$quantiles[,5]
rownames(su.out) <- rownames(summary(dlp.out)$quantiles)
xtable(su.out, digits=3)

# Create vectors of posterior means.
income.vec <- su.out[2:11,1]
lpop.vec <- su.out[12:21,1]
ethfrac.vec <- su.out[22:30,1]
instab.vec <- su.out[31:40,1]
lag.val <- 1:10
# Create matrices for the 95% credible interval.  
income.mat <- su.out[2:11,3:4]
income.mat <- cbind(lag.val, income.mat[,1], lag.val, income.mat[,2])
lpop.mat <- su.out[12:21,3:4]
lpop.mat <- cbind(lag.val, lpop.mat[,1], lag.val, lpop.mat[,2])
lag.val <- 1:9
ethfrac.mat <- su.out[22:30,3:4]
ethfrac.mat <- cbind(lag.val, ethfrac.mat[,1], lag.val, ethfrac.mat[,2])
lag.val <- 1:10
instab.mat <- su.out[31:40,3:4]
instab.mat <- cbind(lag.val, instab.mat[,1], lag.val, instab.mat[,2])
lag.val <- c(1:10)
# Plot all four in one figure.
par(mfrow=c(2,2))
# Income.
par(mar=c(5,5,1,1))
plot(lag.val, income.vec, type="b", pch=19, lwd=3, lty=3, xlab="Lag Value", 
   ylab="Income Posterior Means", ylim=c(-1,1))
abline(h=0, lwd=2, col="grey")
for(i in 1:nrow(income.mat)){
  segments(income.mat[i,1],income.mat[i,2],income.mat[i,3],income.mat[i,4], 
     lwd=2) 
}
# Population.
par(mar=c(5,5,1,1))
plot(lag.val, lpop.vec, type="b", pch=19, lwd=3, lty=3, xlab="Lag Value", 
   ylab="Log Population Posterior Means", ylim=c(min(lpop.mat[,2]),
   max(lpop.mat[,4])))
abline(h=0, lwd=2, col="grey")
for(i in 1:nrow(lpop.mat)){
  segments(lpop.mat[i,1],lpop.mat[i,2],lpop.mat[i,3],lpop.mat[i,4], lwd=2) 
}
# Ethnic fractionalization. 
par(mar=c(5,5,1,1))
lag.val <- c(1:9)
plot(lag.val, ethfrac.vec, type="b", pch=19, lwd=3, lty=3, xlab="Lag Value", 
   ylab="Eth. Fract. Posterior Means", ylim=c(min(ethfrac.mat[,2]), 
   max(ethfrac.mat[,4])))
abline(h=0, lwd=2, col="grey")
for(i in 1:nrow(ethfrac.mat)){
  segments(ethfrac.mat[i,1],ethfrac.mat[i,2],ethfrac.mat[i,3],ethfrac.mat[i,4], 
     lwd=2) 
}
# Political instability.
lag.val <- 1:10
par(mar=c(5,5,1,1))
plot(lag.val, instab.vec, type="b", pch=19, lwd=3, lty=3, xlab="Lag Value", 
   ylab="Instability  Posterior Means", ylim=c(min(instab.mat[,2]), 
   max(instab.mat[,4])))
abline(h=0, lwd=2, col="grey")
for(i in 1:nrow(instab.mat)){
  segments(instab.mat[i,1], instab.mat[i,2], instab.mat[i,3], instab.mat[i,4], 
     lwd=2) 
}

## Table 6: Modernization atheoretic lags model
formula <- demd ~ income + lpop + ethfrac + relfrac + instab + demd_L1 + oil + ncontig     
polity.out <- alreg(formula=formula,
			data=data, model="probit", method="pooled", lags=c(1,2,3,4,5), 
			time.var = "year", by="ccode", min.lag = 1, max.lag = 15,
			ml.type ="Chib95", b0 = 0, B0 = .01,
			mod.prob.method = "BF")

## Figure 5: Compare Model Results Across Applications			
# Reshape the distributed lag modernization example into matrices for plotting.
income.mat[,1] <- income.vec
lpop.mat[,1] <- lpop.vec
ethfrac.mat[,1] <- ethfrac.vec
instab.mat[,1] <- instab.vec
income.mat <- income.mat[,-3]
lpop.mat <- lpop.mat[,-3]
ethfrac.mat <- ethfrac.mat[,-3]
instab.mat <- instab.mat[,-3]

# Generate vectors of baseline and BMA posterior means and credible intervals 
# for plotting, then store all three in a matrix for plotting.
income.rf <- c(summary(rfp.out)$statistics[3,1], summary(rfp.out)$quantiles[3,1], 
   summary(rfp.out)$quantiles[3,5])
income.bma <- c(polity.out$bma[2,1], polity.out$bma[2,3], polity.out$bma[2,4])
income.mat <- rbind(income.rf, income.mat, income.bma)
lpop.rf <- c(summary(rfp.out)$statistics[4,1], summary(rfp.out)$quantiles[4,1], 
   summary(rfp.out)$quantiles[4,5])
lpop.bma <- c(polity.out$bma[3,1], polity.out$bma[3,3], polity.out$bma[3,4])
lpop.mat <- rbind(lpop.rf, lpop.mat, lpop.bma)
ethfrac.rf <- c(summary(rfp.out)$statistics[5,1], summary(rfp.out)$quantiles[5,1], 
   summary(rfp.out)$quantiles[5,5])
ethfrac.bma <- c(polity.out$bma[4,1], polity.out$bma[4,3], polity.out$bma[4,4])
ethfrac.mat <- rbind(ethfrac.rf, ethfrac.mat, ethfrac.bma)
instab.rf <- c(summary(rfp.out)$statistics[7,1], summary(rfp.out)$quantiles[7,1], 
   summary(rfp.out)$quantiles[7,5])
instab.bma <- c(polity.out$bma[6,1], polity.out$bma[6,3], polity.out$bma[6,4])
instab.mat <- rbind(instab.rf, instab.mat, instab.bma)

# Plot Income.
par(mfrow=c(1,1))
plot(c(1:(nrow(income.mat))), income.mat[,1], ylim=c(min(income.mat), 
   max(income.mat)), pch=19, xaxt="n", ylab="Posterior Mean", 
   xlab="Measure", main="Income")
axis(1, c(1:(nrow(income.mat))), labels=c("RF", "DL t-1", "DL t-2", "DL t-3",
   "DL t-4", "DL t-5", "DL t-6", "DL t-7", "DL t-8", "DL t-9", "DL t-10",  "BMA"))
segments(c(1:(nrow(income.mat))), income.mat[,2], c(1:(nrow(income.mat))), 
   income.mat[,3], lwd=2)
abline(h=0, lty=2)     
abline(v=1.5)     
abline(v=11.5)

# Plot Population.
par(mfrow=c(1,1))
plot(c(1:(nrow(lpop.mat))), lpop.mat[,1], ylim=c(min(lpop.mat), max(lpop.mat)), 
   pch=19, xaxt="n", ylab="Posterior Mean", xlab="Measure", 
   main="Population (Logged)")
axis(1, c(1:(nrow(lpop.mat))), labels=c("RF", "DL t-1", "DL t-2", "DL t-3", 
   "DL t-4", "DL t-5", "DL t-6", "DL t-7", "DL t-8", "DL t-9", "DL t-10",  "BMA"))
segments(c(1:(nrow(lpop.mat))), lpop.mat[,2], c(1:(nrow(lpop.mat))), 
   lpop.mat[,3], lwd=2)
abline(h=0, lty=2)     
abline(v=1.5)     
abline(v=11.5)

# Plot Ethnic Fractionalization.
par(mfrow=c(1,1))
plot(c(1:(nrow(ethfrac.mat))), ethfrac.mat[,1], ylim=c(min(ethfrac.mat), 
   max(ethfrac.mat)), pch=19, xaxt="n", ylab="Posterior Mean", xlab="Measure", 
   main = "Ethnic Fract.")
axis(1, c(1:(nrow(ethfrac.mat))), labels=c("RF", "DL t-1", "DL t-2", "DL t-3", 
   "DL t-4", "DL t-5", "DL t-6", "DL t-7", "DL t-8", "DL t-9",  "BMA"))
segments(c(1:(nrow(ethfrac.mat))), ethfrac.mat[,2], c(1:(nrow(ethfrac.mat))), 
   ethfrac.mat[,3], lwd=2)
abline(h=0, lty=2)     
abline(v=1.5)     
abline(v=10.5)

# Plot Political Instability.
par(mfrow=c(1,1))
plot(c(1:(nrow(instab.mat))), instab.mat[,1], ylim=c(min(instab.mat), 
   max(instab.mat)), pch=19, xaxt="n", ylab="Posterior Mean", xlab="Measure",
   main="Instability")
axis(1, c(1:(nrow(instab.mat))), labels=c("RF", "DL t-1", "DL t-2", "DL t-3", 
   "DL t-4", "DL t-5", "DL t-6", "DL t-7", "DL t-8", "DL t-9", "DL t-10",  "BMA"))
segments(c(1:(nrow(instab.mat))), instab.mat[,2], c(1:(nrow(instab.mat))), 
   instab.mat[,3], lwd=2)
abline(h=0, lty=2)     
abline(v=1.5)     
abline(v=11.5)

# Then plot the posterior model probabilities     
plot(polity.out$postprob, type="b", pch=4, lwd=2, xlab = "Feasible Set", 
   ylab = "Posterior Model Probability")

#sink("log.txt")

